Same as Dataset 2f but sample size increased to 100. All different
interaction matrices
1. Load Library
library(tidyverse)
library(igraph)
library(NBZIMM)
library(SpiecEasi)
library(LIMON)
library(here)
library(lme4)
library(Matrix)
library(tscount)
library(patchwork)
library(MASS)
library(matrixcalc)
library(gridExtra)
library(devtools)
library(miaSim)
library(reshape2)
library(corpcor)
2. Simulate Counts
Simulate a different starting covariance matrix and counts per
subject. ie each one has their own seed
# Initalize an empty list
####################################################
count_tables1 <- list()
# The Loop
####################################################
# Set the indices
for (i in 1:100) {
# Set the seed for each subject
set.seed(12345 + i)
# Create Starting covariance matrix
n <- 50^2
lower_values <- runif(n / 2, min = -0.8, max = -0.5)
upper_values <- runif(n / 2, min = 0.5, max = 0.8)
user_interactions <- c(lower_values, upper_values)
user_interactions <- sample(user_interactions)
user_A <- randomA(
n_species = 50,
interactions = user_interactions,
diagonal = -1.0,
connectance = 0.5,
scale_off_diagonal = 0.2,
symmetric = TRUE)
# Heatmap of the covariance matrix
color_palette <- colorRampPalette(c("red", "white", "blue"))(20)
heatmap(
user_A,
Colv = NA,
Rowv = NA,
col = color_palette,
zlim = c(-1, 1),
main = paste("Subject", i)
)
# Run the simulation
tse_glv <- simulateGLV(n_species = 50,
A = user_A,
t_start = 0,
t_store = 4,
stochastic = FALSE,
norm = FALSE,
error_variance = 0.1)
# Get the count table
sim_data <- tse_glv@assays@data@listData[["counts"]]
# Store the count table in the list
count_tables1[[i]] <- t(sim_data)
}




































































































# Merge together
####################################################
# Combine the count tables
combined_count_table <- do.call(rbind, count_tables1)
# Add Rownames
rownames(combined_count_table) <- paste0("Sbj", rep(1:100,
each = nrow(count_tables1[[1]])), "_Time", 1:4)
3. Create fake metadata
Make Metadata and merge with the count data
# Set seed
set.seed(12345)
# Df 1 is Metadata
########################################################
meta_data <- expand.grid(Time = 1:4,ID = 1:100)
rownames(meta_data) <- rownames(combined_count_table)
# Df 2 is Metadata merged with Counts
########################################################
#Round off and increase
combined_count_table <- as.data.frame(combined_count_table + abs(min(combined_count_table)))
combined_count_table <- (combined_count_table)*10
meta_counts <- base::merge(meta_data, combined_count_table, by ="row.names", all = TRUE)
meta_counts <- column_to_rownames(meta_counts, "Row.names")
Plot the data to look at the differences among subjects
Graphs to Check
# Individual Species Plots
########################################################
# Pivot to long data
count_long <- tidyr::pivot_longer(meta_counts, cols = starts_with("sp"), names_to = "Species")
# Plot the data
count_long %>%
ggplot(aes(x = Time, y = value, colour = as.factor(ID),
group = as.factor(ID), linetype = as.factor(ID))) +
geom_line() +
geom_point() +
geom_jitter() +
ylab("Count") +
labs(linetype = "ID", color = "ID") +
facet_wrap(~ Species) + # Create a panel for each species
theme(legend.position = "none") +
ggtitle("Time Series Data")

# Plot only two
########################################################
# Pivot to long data
metacounts_filt <- meta_counts[,c(1:6)]
count_long <- tidyr::pivot_longer(metacounts_filt, cols = starts_with("sp"), names_to = "Species")
# Plot the data
count_long %>%
ggplot(aes(x = Time, y = value, colour = as.factor(ID),
group = as.factor(ID), linetype = as.factor(ID))) +
geom_line() +
geom_point() +
geom_jitter() +
ylab("Count") +
labs(linetype = "ID", color = "ID") +
facet_wrap(~ Species) + # Create a panel for each species
theme_linedraw() +
theme(axis.text.x = element_text(family = "arial",color = "black", size = 14),
axis.text.y = element_text(family = "arial",color = "black", size = 14),
axis.title.x = element_text(family = "arial",color = "black", size = 14),
axis.title.y = element_text(family = "arial",color = "black", size = 14),
legend.position = "none",
strip.text = element_text(face = "bold", family = "arial", color = "white", size = 15))

# Distribution of counts
########################################################
hist(as.matrix(combined_count_table), breaks = 100, main = "Distribution of GLV Data", xlab = "Counts")

4. Save the data
Save the Counts
# Make more like counts
meta_counts[,3:52] <- round(meta_counts[,3:52])
write.csv(meta_counts, here("Data","GLV_SimData", "Dataset_2f", "GLV_100.csv"))
# Save just timepoint 2
meta_counts_filt <- meta_counts %>% filter(Time == 2)
write.csv(meta_counts_filt, here("Data","GLV_SimData", "Dataset_2f", "GLV_100_T2.csv"))
LS0tCnRpdGxlOiAiRGF0YXNldCAyZiAxMDAgU2ltdWxhdGlvbiIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKU2FtZSBhcyBEYXRhc2V0IDJmIGJ1dCBzYW1wbGUgc2l6ZSBpbmNyZWFzZWQgdG8gMTAwLiBBbGwgZGlmZmVyZW50IGludGVyYWN0aW9uIG1hdHJpY2VzCgojIDEuIExvYWQgTGlicmFyeQoqKioKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGlncmFwaCkKbGlicmFyeShOQlpJTU0pCmxpYnJhcnkoU3BpZWNFYXNpKQpsaWJyYXJ5KExJTU9OKQpsaWJyYXJ5KGhlcmUpCmxpYnJhcnkobG1lNCkKbGlicmFyeShNYXRyaXgpCmxpYnJhcnkodHNjb3VudCkKbGlicmFyeShwYXRjaHdvcmspCmxpYnJhcnkoTUFTUykKbGlicmFyeShtYXRyaXhjYWxjKQpsaWJyYXJ5KGdyaWRFeHRyYSkKbGlicmFyeShkZXZ0b29scykKbGlicmFyeShtaWFTaW0pCmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkoY29ycGNvcikKYGBgCgoKCiMgMi4gU2ltdWxhdGUgQ291bnRzCgpTaW11bGF0ZSBhIGRpZmZlcmVudCBzdGFydGluZyBjb3ZhcmlhbmNlIG1hdHJpeCBhbmQgY291bnRzIHBlciBzdWJqZWN0LiBpZSBlYWNoIG9uZSBoYXMgdGhlaXIgb3duIHNlZWQKYGBge3J9CiMgSW5pdGFsaXplIGFuIGVtcHR5IGxpc3QKIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwpjb3VudF90YWJsZXMxIDwtIGxpc3QoKQoKIyBUaGUgTG9vcAojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCgojIFNldCB0aGUgaW5kaWNlcwpmb3IgKGkgaW4gMToxMDApIHsKICAKICAjIFNldCB0aGUgc2VlZCBmb3IgZWFjaCBzdWJqZWN0CiAgc2V0LnNlZWQoMTIzNDUgKyBpKSAgCiAgCiAgIyBDcmVhdGUgU3RhcnRpbmcgY292YXJpYW5jZSBtYXRyaXgKICBuIDwtIDUwXjIgIAogIGxvd2VyX3ZhbHVlcyA8LSBydW5pZihuIC8gMiwgbWluID0gLTAuOCwgbWF4ID0gLTAuNSkgIAogIHVwcGVyX3ZhbHVlcyA8LSBydW5pZihuIC8gMiwgbWluID0gMC41LCBtYXggPSAwLjgpCiAgdXNlcl9pbnRlcmFjdGlvbnMgPC0gYyhsb3dlcl92YWx1ZXMsIHVwcGVyX3ZhbHVlcykKICB1c2VyX2ludGVyYWN0aW9ucyA8LSBzYW1wbGUodXNlcl9pbnRlcmFjdGlvbnMpCiAgdXNlcl9BICA8LSByYW5kb21BKAogICAgbl9zcGVjaWVzID0gNTAsCiAgICBpbnRlcmFjdGlvbnMgPSB1c2VyX2ludGVyYWN0aW9ucywKICAgIGRpYWdvbmFsID0gLTEuMCwKICAgIGNvbm5lY3RhbmNlID0gMC41LAogICAgc2NhbGVfb2ZmX2RpYWdvbmFsID0gMC4yLAogICAgc3ltbWV0cmljID0gVFJVRSkKICAKICAjIEhlYXRtYXAgb2YgdGhlIGNvdmFyaWFuY2UgbWF0cml4CiAgY29sb3JfcGFsZXR0ZSA8LSBjb2xvclJhbXBQYWxldHRlKGMoInJlZCIsICJ3aGl0ZSIsICJibHVlIikpKDIwKQogIGhlYXRtYXAoCiAgICB1c2VyX0EsCiAgICBDb2x2ID0gTkEsCiAgICBSb3d2ID0gTkEsCiAgICBjb2wgPSBjb2xvcl9wYWxldHRlLAogICAgemxpbSA9IGMoLTEsIDEpLCAKICAgIG1haW4gPSBwYXN0ZSgiU3ViamVjdCIsIGkpCiAgKQogIAogICMgUnVuIHRoZSBzaW11bGF0aW9uCiAgdHNlX2dsdiA8LSBzaW11bGF0ZUdMVihuX3NwZWNpZXMgPSA1MCwKICAgICAgICAgICAgICAgICAgICAgICAgIEEgPSB1c2VyX0EsCiAgICAgICAgICAgICAgICAgICAgICAgICB0X3N0YXJ0ID0gMCwgCiAgICAgICAgICAgICAgICAgICAgICAgICB0X3N0b3JlID0gNCwKICAgICAgICAgICAgICAgICAgICAgICAgIHN0b2NoYXN0aWMgPSBGQUxTRSwKICAgICAgICAgICAgICAgICAgICAgICAgIG5vcm0gPSBGQUxTRSwKICAgICAgICAgICAgICAgICAgICAgICAgIGVycm9yX3ZhcmlhbmNlID0gMC4xKQoKICAjIEdldCB0aGUgY291bnQgdGFibGUKICBzaW1fZGF0YSA8LSB0c2VfZ2x2QGFzc2F5c0BkYXRhQGxpc3REYXRhW1siY291bnRzIl1dCiAgCiAgIyBTdG9yZSB0aGUgY291bnQgdGFibGUgaW4gdGhlIGxpc3QKICBjb3VudF90YWJsZXMxW1tpXV0gPC0gdChzaW1fZGF0YSkKfQoKIyBNZXJnZSB0b2dldGhlcgojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCgojIENvbWJpbmUgdGhlIGNvdW50IHRhYmxlcwpjb21iaW5lZF9jb3VudF90YWJsZSA8LSBkby5jYWxsKHJiaW5kLCBjb3VudF90YWJsZXMxKQoKIyBBZGQgUm93bmFtZXMKcm93bmFtZXMoY29tYmluZWRfY291bnRfdGFibGUpIDwtIHBhc3RlMCgiU2JqIiwgcmVwKDE6MTAwLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGVhY2ggPSBucm93KGNvdW50X3RhYmxlczFbWzFdXSkpLCAiX1RpbWUiLCAxOjQpCmBgYAoKCgoKIyAzLiBDcmVhdGUgZmFrZSBtZXRhZGF0YSAgCgoKTWFrZSBNZXRhZGF0YSBhbmQgbWVyZ2Ugd2l0aCB0aGUgY291bnQgZGF0YQpgYGB7cn0KIyBTZXQgc2VlZApzZXQuc2VlZCgxMjM0NSkKIyBEZiAxIGlzIE1ldGFkYXRhCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCm1ldGFfZGF0YSA8LSAgZXhwYW5kLmdyaWQoVGltZSA9IDE6NCxJRCA9IDE6MTAwKQpyb3duYW1lcyhtZXRhX2RhdGEpIDwtIHJvd25hbWVzKGNvbWJpbmVkX2NvdW50X3RhYmxlKQoKCgojIERmIDIgaXMgTWV0YWRhdGEgbWVyZ2VkIHdpdGggQ291bnRzCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCiNSb3VuZCBvZmYgYW5kIGluY3JlYXNlCmNvbWJpbmVkX2NvdW50X3RhYmxlIDwtIGFzLmRhdGEuZnJhbWUoY29tYmluZWRfY291bnRfdGFibGUgKyBhYnMobWluKGNvbWJpbmVkX2NvdW50X3RhYmxlKSkpCmNvbWJpbmVkX2NvdW50X3RhYmxlIDwtIChjb21iaW5lZF9jb3VudF90YWJsZSkqMTAKbWV0YV9jb3VudHMgPC0gYmFzZTo6bWVyZ2UobWV0YV9kYXRhLCBjb21iaW5lZF9jb3VudF90YWJsZSwgYnkgPSJyb3cubmFtZXMiLCBhbGwgPSBUUlVFKQptZXRhX2NvdW50cyA8LSBjb2x1bW5fdG9fcm93bmFtZXMobWV0YV9jb3VudHMsICJSb3cubmFtZXMiKQpgYGAKCgpQbG90IHRoZSBkYXRhIHRvIGxvb2sgYXQgdGhlIGRpZmZlcmVuY2VzIGFtb25nIHN1YmplY3RzCgpHcmFwaHMgdG8gQ2hlY2sKYGBge3J9CiMgSW5kaXZpZHVhbCBTcGVjaWVzIFBsb3RzCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCiMgUGl2b3QgdG8gbG9uZyBkYXRhCmNvdW50X2xvbmcgPC0gdGlkeXI6OnBpdm90X2xvbmdlcihtZXRhX2NvdW50cywgY29scyA9IHN0YXJ0c193aXRoKCJzcCIpLCBuYW1lc190byA9ICJTcGVjaWVzIikKCiMgUGxvdCB0aGUgZGF0YQpjb3VudF9sb25nICU+JQogIGdncGxvdChhZXMoeCA9IFRpbWUsIHkgPSB2YWx1ZSwgY29sb3VyID0gYXMuZmFjdG9yKElEKSwKICAgICAgICAgICAgIGdyb3VwID0gYXMuZmFjdG9yKElEKSwgbGluZXR5cGUgPSBhcy5mYWN0b3IoSUQpKSkgKwogIGdlb21fbGluZSgpICsgCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2ppdHRlcigpICsKICB5bGFiKCJDb3VudCIpICsKICBsYWJzKGxpbmV0eXBlID0gIklEIiwgY29sb3IgPSAiSUQiKSArCiAgZmFjZXRfd3JhcCh+IFNwZWNpZXMpICsgICMgQ3JlYXRlIGEgcGFuZWwgZm9yIGVhY2ggc3BlY2llcwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKwogIGdndGl0bGUoIlRpbWUgU2VyaWVzIERhdGEiKQoKCiMgUGxvdCBvbmx5IHR3bwojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwojIFBpdm90IHRvIGxvbmcgZGF0YQptZXRhY291bnRzX2ZpbHQgPC0gbWV0YV9jb3VudHNbLGMoMTo2KV0KY291bnRfbG9uZyA8LSB0aWR5cjo6cGl2b3RfbG9uZ2VyKG1ldGFjb3VudHNfZmlsdCwgY29scyA9IHN0YXJ0c193aXRoKCJzcCIpLCBuYW1lc190byA9ICJTcGVjaWVzIikKCiMgUGxvdCB0aGUgZGF0YQpjb3VudF9sb25nICU+JQogIGdncGxvdChhZXMoeCA9IFRpbWUsIHkgPSB2YWx1ZSwgY29sb3VyID0gYXMuZmFjdG9yKElEKSwKICAgICAgICAgICAgIGdyb3VwID0gYXMuZmFjdG9yKElEKSwgbGluZXR5cGUgPSBhcy5mYWN0b3IoSUQpKSkgKwogIGdlb21fbGluZSgpICsgCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2ppdHRlcigpICsKICB5bGFiKCJDb3VudCIpICsKICBsYWJzKGxpbmV0eXBlID0gIklEIiwgY29sb3IgPSAiSUQiKSArCiAgZmFjZXRfd3JhcCh+IFNwZWNpZXMpICsgICMgQ3JlYXRlIGEgcGFuZWwgZm9yIGVhY2ggc3BlY2llcwogIHRoZW1lX2xpbmVkcmF3KCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGZhbWlseSA9ICJhcmlhbCIsY29sb3IgPSAiYmxhY2siLCBzaXplID0gMTQpLAogICAgICAgIGF4aXMudGV4dC55ID0gZWxlbWVudF90ZXh0KGZhbWlseSA9ICJhcmlhbCIsY29sb3IgPSAiYmxhY2siLCBzaXplID0gMTQpLAogICAgICAgIGF4aXMudGl0bGUueCA9IGVsZW1lbnRfdGV4dChmYW1pbHkgPSAiYXJpYWwiLGNvbG9yID0gImJsYWNrIiwgc2l6ZSA9IDE0KSwKICAgICAgICBheGlzLnRpdGxlLnkgPSBlbGVtZW50X3RleHQoZmFtaWx5ID0gImFyaWFsIixjb2xvciA9ICJibGFjayIsIHNpemUgPSAxNCksCiAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiLAogICAgICAgIHN0cmlwLnRleHQgPSBlbGVtZW50X3RleHQoZmFjZSA9ICJib2xkIiwgZmFtaWx5ID0gImFyaWFsIiwgY29sb3IgPSAid2hpdGUiLCBzaXplID0gMTUpKQoKCiMgRGlzdHJpYnV0aW9uIG9mIGNvdW50cwojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwpoaXN0KGFzLm1hdHJpeChjb21iaW5lZF9jb3VudF90YWJsZSksIGJyZWFrcyA9IDEwMCwgbWFpbiA9ICJEaXN0cmlidXRpb24gb2YgR0xWIERhdGEiLCB4bGFiID0gIkNvdW50cyIpCmBgYAoKCiMgNC4gU2F2ZSB0aGUgZGF0YQoKU2F2ZSB0aGUgQ291bnRzCmBgYHtyfQojIE1ha2UgbW9yZSBsaWtlIGNvdW50cwptZXRhX2NvdW50c1ssMzo1Ml0gPC0gcm91bmQobWV0YV9jb3VudHNbLDM6NTJdKQp3cml0ZS5jc3YobWV0YV9jb3VudHMsIGhlcmUoIkRhdGEiLCJHTFZfU2ltRGF0YSIsICJEYXRhc2V0XzJmIiwgIkdMVl8xMDAuY3N2IikpCgojIFNhdmUganVzdCB0aW1lcG9pbnQgMgptZXRhX2NvdW50c19maWx0IDwtIG1ldGFfY291bnRzICU+JSBmaWx0ZXIoVGltZSA9PSAyKQp3cml0ZS5jc3YobWV0YV9jb3VudHNfZmlsdCwgaGVyZSgiRGF0YSIsIkdMVl9TaW1EYXRhIiwgIkRhdGFzZXRfMmYiLCAiR0xWXzEwMF9UMi5jc3YiKSkKYGBgCgoKCgo=